Random Diffusion Model with Structure Corrections 
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Tlie random diffusion model is a continuum model for a conserved scalar density field driven 
by diffusive dynamics where the bare diffusion coefficient is density dependent. We generalize the 
model from one with a sharp wavenumber cutoff to one with a more natural large-wavenumber 
cutoff. We investigate whether the features seen previously - namely a slowing down of the system 
and the development of a prepeak in the dynamic structure factor at a wavenumber below the first 
structure peak - survive in this model. A method for extracting information about a hidden prepeak 
in experimental data is presented. 
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I. INTRODUCTION 

In this work we study a generalization of the random 
diffusion model (RDM) introduced by Mazenko in Ref. [ll, 
henceforth referred to as RDMI. The RDM is a model for 
nonlinear diffusion in colloidal systems. It is formulated 
in terms of a conserved density field and, in its simplest 
form, can be chosen to have Gaussian static statistics. 
The model can be motivated as a continuum general- 
ization of the facilitated spin models of glassy dynamics 
[1^ 8] . These models have a density dependent kinetic co- 
efficient which leads to slowing down in dense systems. 
The RDM similarly has a density dependent bare diffu- 
sion coefficient. In RDMI, this bare diffusion coefficient 
led to a slowing down as the density was increased. The 
more realistic RDM presented here exhibits similar be- 
havior. 

There has been much speculation but relatively few 
solid results in establishing the existence of a mode- 
coupling theory (MOT) ergodic-nonergodic (ENE) tran- 
sition in the field-theoretic models of the liquid-glass 
transition. The RDM is a candidate for the simplest such 
model, however this point will not be a focus of this pa- 
per. In RDMI, the model was shown to undergo an ENE 
transition at one-loop order, but not at two-loop order. 
The RDM is part of a larger class of models with den- 
sity dependent bare diffusion coefficients including those 
discussed by Dean Q , Kawasaki and Miyazima [lO] , and 
Miyasaki and Reichman fTT'|. 

In RDMI, the Fokker-Planck formalism was used to 
generate a self-consistent perturbation expansion for the 
memory function. This memory function could be in- 
serted into the kinetic equation which describes the time- 
evolution of the dynamic correlation function, the phys- 
ical observable of interest. In the simplest realization of 
this model, one has a single control parameter, g, which 
controls the density dependence. RDMI treated a course- 
grained system in what was termed the structureless ap- 
proximation; the static structure factor of the system was 
assumed constant up to a wavenumber cutoff A. 

Solving the kinetic equation revealed the following: 



1. As 5 increases, the system slows down. 

2. A new peak, termed the prepeak, develops at a 
wavenumber qq away from zero. The form of this 
peak can be fit to a Gaussian 



(1) 



3. 



where the amplitude A decreases with time and the 
peak width narrows with time. This form 

reveals a new growing length in the problem (v^) 
and fits of the the amplitude change smoothly from 
an exponential form at zero coupling to a power law 
form near a critical value. 

Above this critical value, the system becomes un- 
stable and the prepeak amplitude grows without 
bound. Such behavior is clearly unphysical and as 
such, the model breaks down. 



In this paper, we generalize the random diffusion model 
to a more realistic form where wavenumber integrations 
are cut off naturally in the theory. This requires incor- 
porating a more realistic static structure factor into the 
theory. Here, the integrals are naturally cut off by the 
large wavenumber decay of the direct correlation func- 
tion. Our new model depends on two parameters; we 
have a coupling constant analogous to that seen in the 
structureless case, but our model also depends explicitly 
on density. 

For this new model we find the following: 

1. As the coupling and density increase, the system 
again slows down. 

2. A prepeak forms, now located between q = and 
the first structure peak. This peak and the peaks 
due to the static structure factor can be fit to Gaus- 
sians with decaying amplitude and narrowing width 
just as before and the amplitude again changes 
from an exponential to a power law form as the 
system approaches the critical crossover. 
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3. Above the critical values of density and coupling 
constant, we again find that the prepeak grows 
without bound. Whereas previously this separation 
between stable and unstable growth was character- 
ized by a single number, our model has a separation 
given by a critical line in coupling constant-density 
parameter space. 

This improved model therefore preserves the features 
seen in RDMI. 

In this work, we again discuss in detail the slowing 
down of the system and the interesting feature of the 
dynamically generated prepeak. Though the develop- 
ment of this prepeak in both RDMI and this paper is 
surprising, it is worth pursuing. The random diffusion 
model is simple, yet incorporates nonlinearities known to 
be present and the perturbation is straightforward and 
without the "tricks" sometimes used to force things into 
mode-coupling form. It is dynamically generated instead 
of arising via the static structure factor and we find here 
that it is robust in the choice of said statics; what one 
may have worried was an artifact of a coarse approxi- 
mation in RDMI is shown here to survive in the same 
form. 

The final feature of the model, the unstable growth of 
the system above critical values of the model parameters, 
is a feature which reveals the breakdown of the model. 
In this region, the system likely wants to nucleate, but 
no terms are present in the model to provide this stabi- 
lization. We discuss some loose features of the instability 
and then propose the terms, motivated by density func- 
tion theory, which may cutoff the growth. 

Because our model uses a realistic structure factor, we 
can draw a closer connection to experimental results than 
was possible in RDMI. In the last section of the paper we 
present a method by which one can extract information 
about a possible prepeak hidden in experimental mea- 
surements of the dynamic structure factor. Though few, 
some experimental allusions to prepeaks exist and com- 
parison is desired. 



II. THE RANDOM DIFFUSION MODEL 

A. Introduction 

Let us introduce the random diffusion model in the 
Fokker-Planck context. For a fundamental density field 
0(x), wc take an effective Hamiltonian quadratic in <p 
given by 



(2) 



where ^0(x) = 0(x) — n and n — ((/>). We will study 
density time correlations through the equilibrium inter- 
mediate dynamic structure factor given by 

^^(qi,q2;t) - (0(q2,t)0(qi,O)) = (0(q2)e-^^V(qi)) 
= (27r)'^5(qi+q2)C(qi;t) (3) 



where <^(q) is the Fourier transform of the fundamental 
field ^0, averages are given by = / 'D{(j))W^f{(j)) 

with an equilibrium probability distribution of the usual 
form, = e~^'^'t'/Z, and is the adjoint Fokker- 
Planck operator defined below. We will also make fre- 
quent use of the static (equal-time) correlation function 
given by 

C(qi,q2) = C(qi,q2;i - 0) = (2^)'^^(qi + q2)C'(qi)(4) 

which is related to the static susceptibility appearing in 
Eq. @ by 

^(qi)=rV(qi)- (5) 
The adjoint Fokker-Planck operator takes the form 



5Hi 



hT 



^0(x) 



r0(x,y) 



(6) 



where is a transport matrix which incorporates the 
bare diffusion coefficient [11] . In RDMI, F^ was given by 

r^(x,y) = V,V,[i^(0),5(x-y)] (7) 

where the bare diffusion coefficient D{(p) was 

D{(j>) ^ Do + Di<t){x) = D + Di5(j>{x.) (8) 

with D — {D{(f))) — Dq + Din. For this work, however, 
we choose a more general form 

r^(x,y) = V^Vy J d''z(^foix~z)Dfo{y-z) 

+ hix~z)DiScPiz)h{y^z) + ..y (9) 

By relaxing the constraining delta-function and introduc- 
ing the general functions /o and /i, we will have consid- 
erably more freedom in regulating integrals which result 
from application of perturbation theory. For now, we 
leave the functions undetermined, but will later set them 
in such a way as to give the short-time sum rules correctly 
and give the large- wavenumber dependence of the vertex 
in the memory function as in mode-coupling theory. 

Our choice for the transport matrix (both the form 
of RDMI and the more general form here) is motivated 
as a course grained alternative to the microscopic form 
D{(j)) = Do most often used. Density dependent terms 
are expected to play a role in the dynamics and as such, 
have been incorporated here. 

Our model therefore depends on the functions fo, fi 
and X s-nd the constants Do and Di. The effects of ex- 
tending the transport matrix to higher order through ad- 
ditional constants (D2, etc.) and functions (/2, etc.) will 
be addressed in a future paper. 

The physical diffusion coefficient is given by 



Dp = lim D{q) 



(10) 
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where D(q) comes from the Fourier transform of the av- 
erage of the transport matrix: 

D{q)q^ - j d'^xie"l (''^-^=)(r^(xi - X2)). (11) 

Thus, we have 

Dp = Df^{Q) (12) 

where we have used the fact that averages over odd pow- 
ers of 5<j) vanish. 



where F(q, k) is the vertex function. We define our ver- 
tex in light of our new choice for as 

t/(q,k) = ^i?i[/i(<z)q-/i(fc)kx-^(fc) 

+ /i(<z)q • /i(q - k)(q - k)x-i(q - k)](19) 

which has a simple dependence on the function /i. We 
choose this form to mimic the traditional mode-coupling 
theor y fo rm. (For more on traditional MCT, see, e.g. 
Refs. |15| and|l6|.) Now that we have defined this vertex, 
it is set to arbitrarily high order. 



B. Memory Function Formalism 

In order to study the time evolution of the dynamic 
structure factor C(g, t), we organize our theory using the 
memory function technique |14l | . In Laplace transformed 
space, the kinetic equation is given by 

[z + K{q,z)]C{q,z)^C{q) (13) 

where 

C{q,z) = -i dte'''C{q,t) (14) 

and where K{ci,t) = K'''\ci) + K'-'^\ci,t) is the memory 
function decomposed into a static piece and a dynamic 
piece. Inverse Laplace transforming this equation gives 
the form 

^ = ^K^^H<l)C{q,t) 

+ f dsK'^'^\q,t- s)C{q,s). (15) 
Jo 

First, the static memory function is determined (with- 
out approximation) by the equilibrium average 

K(^\q)C{q)=ir'{yM)- (16) 
Using our transport matrix given by Eq. ([S]), we have 

if(^)(q) = ir\^Dfl{q)C-\q) 

= n^Dfl{q)x-\q). (17) 

The dynamic part of the memory function is more com- 
plicated to compute and can be determined via a pertur- 
bation expansion. If we assume that the nonlinearities 
are small, we may expand in powers of Di, the coeffi- 
cient of the density correction to the diffusion coefficient. 
A full derivation is given in RDMI, but at lowest order, 
the dynamic memory function is given by 

i^('^)(q,i)C(<7) 

^-s/ ^[nq,k)FC'(k;i)C(q-k;t) (18) 



C. Structure corrections 

Perturbation expansions inevitably run into wavenum- 
ber integrals such as the one we derived for the dynamic 
memory function, Eq. (|18p . which are potentially diver- 
gent. Care must be taken when forming and evaluat- 
ing the integrals, and one usually must either develop a 
rationale for implementing a large-wavenumber cutoff or 
structure the vertex such that the integrals remain finite. 

Constructing a vertex which enforces convergence can 
be difficult. For example, Ref. [l3| studies the Dean- 
Kawasaki model 0, [l^ applicable to colloids. This work 
analyzes the theory in terms of time-reversal symmetry, 
but this approach ultimately leads to a vertex that is 
not of a standard MCT form and to divergent integrals. 
They propose that the solution may be a cumbersome 
resummation of higher-order diagrams to renormalizc the 
vertex, but leave a full solution to future work. 

For the random diffusion model, divergent integrals 
first appeared in RDMI. The approach in that work was 
to take a sharp cutoff at finite wavenumber A and (for 
most of the paper) to take the static susceptibility to be 
constant such that x^'^il) ~ f- This approach was called 
the structureless approximation because it corresponded 
to a coarse-grained system where the short-distance de- 
grees of freedom (including the first peak in the static 
structure factor) had been integrated out jl3j] . While this 
simplified the problem and kept the integrals finite, one 
has reason to worry whether the results are strongly in- 
fluenced by the nature of the cutoff and the structureless 
nature of the structure factor. It is therefore desirable 
to to impose a different method to regulate the integrals 
which does not raise such concerns. 

To regulate the wavenumber integrals in this work, 
we now define our functions Xj /o, /i introduced ear- 
lier. First, let us restore the structure to the model by 
defining xil) in terms of the full static structure factor, 
S{q), as 

X{q) = npS{q) (20) 

where 
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and Coiq) is the direct correlation function which de- 
cays to zero at large wavenumber. For simplicity, we will 
use the solution to the Percus-Yevick approximation for 
hard-spheres in this paper. (See e.g. Refs. [l9l and I20I for 
details.) In principle, one could take any realistic ana- 
lytic approximation for S{q) or use experimental results. 

Next, the function fo{q) appears in the static mem- 
ory function K'^'^^q) and the physical diffusion coefficient 
Dp{q). Since we have no regularization constraints on ei- 
ther of these functions, let us take foiq) ~ 1 which gives 



iq^Df^(q) iq^D 



x{q) 



C{q)P 



and 



(22) 



(23) 



Finally, the function fi{q) appears only in the ver- 
tex (and therefore in the dynamic memory function 
K^'^\q,t)). As discussed above, we want to choose fi{q) 
so that the memory function integral remains finite as we 
take the cutoff to infinity. Let us reason our way through 
the appropriate choice. 

As it stands, the long wavenumber behavior of the ver- 
tex (beside the explicit dependence) is governed by 
X'^il) which approaches unity at as g — 00. Let us in- 
stead let the vertex go as the direct correlation function 
Cniq) which approaches zero in the same limit. Thus, 
we want fi{q) ~ xi(l)C^D{q)- Normalizing fi{q) to unity 
at g = 0, we therefore have 



x{q)CD{q) ^ S{q)nCD{q) 
x(0)Cz,(0) S{Q)nCD{Qy 



(24) 



Note that if we now put fi{q) back into the bare dif- 
fusion coefficient (Eq. ^) and take the large-g (short 
distance) limit, the density dependent Di term goes to 
zero returning the microscopic result D{(j)) — D. 



This now gives the vertex 



V{q,k) = 



iDi S{q)nCD{q) 



2n/3 (S'(0)nCi)(0))2 
X q-[k7iCD(k) + (q-k)nCz3(q-k)].(25) 

This is the traditional MCT form which decays to zero 
as either g or /c — ^ 00. Such behavior, we will see, is 
sufficient to keep the memory function finite without im- 
posing a finite integral cutoff. 

To simplify the kinetic equation let us now introduce 
dimensionless variables. The dimensionless wavenumber, 
time and density (packing fraction) are given by 



Q = q<^, 



T = Datp^^ 



(26) 



(27) 



and 



11 = (7r/6)ncr^ (28) 
where a is the hard sphere diameter and the dimension- 
less correlation function is given by 



f{Q,T)^C{Q,T)/C{Q). 



(29) 



After inserting the static and kinetic equations in terms 
of dimensionless variables, we have 



dT 



fiQ.T) 



'^(g)677 

-R' I dSN{Q,T-S)f{Q,S) (30) 



where we've defined the simplified memory function 



J 



N{Q,T) = 



1 



i:yS{Q)n-'Cl{Q) 



2 \6r]J S'4(0)n4C4 (0) J (27r)3 
X SiK)f{K,T)S{Q-K)fiQ-K,T) 



Q • KnCoiK) + Q iQ- K)nCi5(Q - K) 



r 



(31) 



and the coupling constant 
Din 



R = 



D 



Din 



Do + nDi 



(32) 



Let us make a few comments on this form. 

First, the kinetic equation now depends only on two 
parameters, the density rj and the coupling constant 
R. Note that the temperature dependence has been 
absorbed into the dimensionless variables. The kinetic 



equation can be solved numerically to get f{Q,T) and 
the general behavior as rj and R are varied can be stud- 
ied. 

Second, the coupling constant R can be either positive 
or negative (as the coefficient Di can be either positive 
or negative). However, since R enters only through R^, 
the overall sign of Di is irrelevant. We choose for our 
purposes the constraint < i? < 1. One can imagine 
combinations involving negative, but large values of Din 
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such that \R\ > 1, but since our perturbation expansion 
assumed small Z^i, we are already overly generous by 
studying R up to unity. 

Finally, note that our form for df{Q,T)/dT diverges 
as as 77 — 0. While this may seem strange, this 
is just a manifestation of the characteristic time scale of 
the system slowing down with density. One can choose to 
rescale the dimensionless time as T — > T' = Tn/dri and 
the explicit divergence disappears. For various reasons, 
we have chosen to use the unsealed time, though if one 
were interested in low density systems, the scaled variable 
would be a more appropriate choice. 
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III. NUMERICAL RESULTS 

A. RDMI: A Review 

Before we solve our model, let us review in more de- 
tail the results from RDMI. In that model, the kinetic 
equation (at lowest order) is given by 



df{Q,T) 
dT 



= -Q'f{Q,T) 

+ Q* dSN{Q,T- S)f{Q,S) (33) 



with 



(2^ 



/(i^,r)/(Q-K,T) (34) 



and 



(35) 



where C is a constant and A is the wavenumber cutoff. 

In the non-interacting case {g = 0), the equation can 
be solved analytically to get the simple form 



/o(Q,T) = e 



(36) 



As the coupling g increases, the decay at large 
wavenumber values slows and a small peak develops. 
This peak (termed the prepeak) can be fit to a Gaus- 
sian of the form 

/p(Q,r)=v4e-^(«3-Q°)' (37) 

where the amplitude A has the time dependence 



A{T) = Ao 



-ET 



(T + TaY 



while B satisfies 



B{T) = Ba{T + TBf 



(38) 



(39) 



where Qo, Aq, Bq, Ta, Tb, a, /3 and E are positive, 
time-independent fit parameters. 



FIG. 1: A plot of f{Q,T) as determined in RDMI for a cou- 
pling g > g* . Each line is equally spaced in time and the 
times are earliest to latest as the peak amplitude goes from 
smallest to largest. Wavenumber space is normalized such 
that the cutoff corresponds with Q=l. f{Q,T)=0 for larger 



As the system approaches the critical coupling g*, the 
peak amplitude goes to a completely algebraic decay with 
time such that — > 0. Above g*, the above analytic 
equations cease to hold; the peak no longer decays at 
long times, but instead grows without bound, rendering 
the system unstable. Examples of the behavior of this 
model can be seen in Fig. [T] 



B. General Results 

For our model, the evolution of the dynamic correlation 
function is given by Eq. ([50]). While similar to the form 
of RDMI, the incorporation of a realistic static structure 
factor leads to a slightly more complicated form. We can 
numerically solve the kinetic equation for a particular set 
of the parameters rj and R. (The details of this solution 
are given in Appendix A.) The upper bound on the value 
of 77 is 0.74 which is the packing fraction of a close-packed 
solid. The upper limit for i? is 1 as discussed at the end 
of the last section. 

Before looking at the complete problem, let us first ex- 
amine the simpler R = case to see some of the features 
which are common to all the solutions of the kinetic equa- 
tion. In this limit there is no dynamic feedback from the 
memory function and the kinetic equation can be solved 
analytically to give 



(40) 



Thus, the correlation function decays exponentially with 
time and its wavenumber dependence is strongly deter- 
mined by the static structure factor S{Q). This is the 
well known de Gennes narrowing form if we take a scaled 
time T' = Ttt/G?]. 
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FIG. 2: Plots of /(Q, T) with 7? = for (a) 77 = 0.3, (b) 
r) = 0.4, (c) ?7 = 0.5 and (d) 77 = 0.6. Times are (from top 
line to bottom line) T^O.OOl, 0.01, and 0.1. 



In Fig. [2] we plot the behavior of f{Q,T) for various 
values of rj. We see that a number of peaks form which 
decay to zero with time. Beside the peak centered at 
Q—0, these peaks are centered at the same wavenumbers 
as the peaks of the static structure factor. Even with 
i? = 0, we see that we have considerable structure in the 
theory. 

We now turn to the case of i? 7^ 0. 

At small 77, say 77 = 0.10 or 0.20, the behavior is nearly 
indistinguishable from the R = case even as R is in- 
creased through all allowed values. As we increase to 
intermediate values of 77, say rj — 0.30, the situation be- 



comes more interesting. In Fig. 3(a) we plot results for 
77 — 0.30 and R = 0.51 and see a slight slowing down of 
the decay. That is, the peaks in the R — 0.51 case persist 
longer before decaying to zero than those in the R — 
case. 

If we stay at these intermediate densities and ratchet 
up to extremely large R, we see an interesting new fea- 
ture; between the Q=0 peak and the first structure peak, 
a new small peak appears. With time, this peak also de- 
cays away to zero and disappears before the first struc- 
ture peak. Because of its position in Q-space, we will call 
this peak the prepeak because it is located at a smaller 
Q value than the first structure peak. As we push the 
coupling up to the limit i? = 1, we find that the prepeak 
persists for longer and longer times, but still decays to 
zero with all the other peaks. (Fig. |3(b)[ ) 

Looking at the behavior of the second and higher struc- 
ture peaks, we find little of interest. These higher order 
peaks mirror the behavior of the first structure peak - 
their decay slows with increasing R - and wc find no new 
features (e.g., higher order prepeaks). This remains true 
for most of the work presented in this paper and we will 
therefore restrict ourselves to discussion of the behav- 
ior at wavenumbers around and below the first structure 
peak unless otherwise noted. 

Moving to 77 = 0.40, we again slowly increase the value 
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FIG. 3: Plots of /(Q,r) for r; = 0.30 with (a) R = 0.51 and 
(b) i? = 1. Times are (from top to bottom) r=0.02, 0.03, 0.05 
and 0.15. Note that a prepeak is developing in the R = 1 plot, 
however both plots decay to zero with time. 



of R. Initially, we clearly see the first structure peak and 
note that it decays to zero with time as in the R — 
case. We plot f(Q,T) for R = 0.52 in Fig. 
ample. Increasing R a bit further to i? = 
we cross over to a new regime. Now, 



4(a; 



as an ex- 



0.56 however, 
the first structure 
peak amplitude initially decreases, then slows down and 
finally begins to increase. (Fig. |4(b)[ ) If we continue 
the solution to long times, we find that this growth has 
no bound and in fact accelerates. The prepeak (which is 
initially only weakly visible) grows faster than the first 
structure peak and the two peaks are of comparable am- 
plitude only after the model and the numerical solution 
break down. We plot the amplitude of the prepeak and 
first structure peak in Fig. [5l including late times where 
the solution has become unphysical. 

We may continue in this manner and look at higher 
values of 77. We see a very general pattern appear. As R 
is increased from zero, there is initially a period where all 
peaks decay to zero with time. The plots look qualita- 
tively like the R = plots, however as R is increased, the 
peaks decay more slowly. We call this region of parame- 
ter space the stable region. Continuing, we pass a critical 
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value of R and see that the amplitudes of the peaks no 
longer decay to zero, but turn around and grow without 
bound. We call this (interesting, yet unphysical) region 
the unstable region and label the critical coupling R* . We 
find that as we increase 77, the value of R* at which we 
cross from the stable to unstable regime decreases. As we 
continue increasing the value of R past R*, we find that 
the peak turns around its growth at earlier and earlier 
times. 

The behavior of the prepeak is quite curious. As we 
increase rj, we find that the prepeak becomes more diffi- 
cult to see in the plots. Recall that at 77 = 0.30, we could 
clearly see the prepeak even though we were in the sta- 
ble region. At high values of 77 we do not see the prepeak 
in the stable region and often do not see the prepeak in 
the unstable region until near the numerical breakdown 
when both the prepeak and first structure amplitudes be- 
come comparable. When both peaks are visible, we note 
that the first structure peak turns around and grows at 
an earlier time than the prepeak, but that the prepeak 
grows at a faster rate. 

Though the prepeak is sometimes not visible from the 
raw data, we can extract some information about the 
prepeak's behavior. We believe that the behavior of the 
prepeak and the first structure peak, while quantitatively 
different, are related and that we do not lose information 
by monitoring just the latter. We discuss this in more 
detail in the sections which follow. 



FIG. 4: Plots of f{Q, T) for r) = 0.40 with (a) R = 0.52 and 
(b) R — 0.56. Times are (from top to bottom at Q = 0) 
r=0.05, 0.1 and 0.3. Note that whereas the first structure 
peak of the R — 0.52 plot decays to zero with time, the peak 
in the R — 0.56 plot grows. 
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FIG. 5: Logarithmic plots of the first structure peak am- 
pfitude (solid line) and prepeak amplitude (dashed line) for 
r] = 0.40, R — 0.56. Note that even though the first structure 
peak amplitude begins to grow first, the prepeak amplitude 
grows at a faster rate and the two are comparable when the 
solution breaks down. 



IV. ANALYSIS 
A. Stable-unstable crossover 

We have seen that our solution has two distinct 
regimes: the stable region and the unstable region. We 
now wish to quantify the crossover from stable to un- 
stable behavior. Because this crossover depends on both 
the coupling constant and the density, we have a critical 
function R*{r]). 

To determine this line, we follow the approach devel- 
oped in RDMI. Let us first hold 77 fixed and vary R. 
If we are in the stable region, the amplitudes of the 
peak heights decrease monotonically with time. In the 
unstable region however, we may identify a time T^m 
where the amplitude of a particular peak reaches its low- 
est point. (See, for example. Fig. [5]) As we approach 
R* from larger R values, we find that the value of Tmm 
approaches infinity as we get arbitrarily close to R* . We 
may therefore find T^m for values of R close to but larger 
than R* and fit these values to a function diverging at 
R*: 



Tn 



{R- R*)^ 



(41) 



If we repeat this technique for a set of ij values, we de- 
termine a corresponding set of R* values along the line 
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FIG. 6: Plot of the critical line R*{rj) separating the stable 
and unstable regions. Each point was determined from a set 
of data using Eq. (|4T|. The solid line is to guide the eye only. 





Unstable 
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R* (rj) which divides the stable and unstable regimes [21| 
Our results for R*{ri) are plotted in Fig. |6l 



With a technique in place, one now may ask which 
peak to use in determining Tmin- We saw in the previ- 
ous section that the prepeak and the first structure peak 
do not reach their minima at the same time and that 
their unstable growth is not equally paced. It is there- 
fore unclear if the results from fitting one set of data will 
match the results from fitting the other. It is also un- 
clear whether there may in fact be a region where we 
have unstable growth in one peak with stable decay in 
the other. A careful analysis, however, reveals that this 
is not the case. Values for R* determined from the pre- 
peak amplitude closely match those determined from first 
structure peak amplitude and we are unable to find any 
region where the stability /instability of one peak does 
not match that of the other. Thus, we chose to use the 
first structure peak amplitude data for determining the 
critical function since it is the most clearly visible in all 
plots. 

We can note several things from the plot in Fig. [HI 
First, we again see that for small to moderate values 
of rj the solution is stable for all studied values of R. 
The critical line does not cross the R — 1 mark until 
some value between rj = 0.30 and rj = 0.35. Second, 
even though a transition to the unstable region at these 
moderate densities seems possible, one should be cautious 
about such a conclusion. Recall that we expect R to be 
a small perturbation parameter. By this reasoning, we 
may wish to restrict ourselves to discussion of such a 
transition at higher densities where R* is indeed small. 
Such higher densities are are those more typical of very 
dense liquids, solids and glasses. 



B. Analytical Fits 

We now see that the dynamic structure factor has the 
general form 

/(Q, T) = /o(g, T) + /pp(g, T) + /„(Q, T) (42) 

where fo{Q,T) is the initial peak centered at Q = 0, 
/pp(Q, T) is the prepeak and fn{Q, T) are the peaks due 
to the structure factor, of which we are only concerned 
in the first {n = 1). Each peak can be modeled as a 
Gaussian such that 



UQ,T) = A,iT)e-^^(^^(Q- 



(43) 



where x may be 0, pp or f depending on the peak in 
question. 

The peak centered at Q = is the simplest to treat. 
For small Q, we may neglect the interaction term and 
we are left with the de Gennes form given earlier by Eq. 



(44) 



Gomparing this to the generic Gaussian form of Eq. (|43l) , 
we have = 1; Qo = and 



Ttt 
6ryS'(0)' 



(45) 



A fit to the data shows excellent agreement. We find, for 
example, that for 77 = 0.55, R = 0.14, Eq. (1451) predicts 
Bo = 102. 38T while a fit yields Bq = 102. 26T. 

Near the prepeak and first structure peak, the situation 
is more complicated. In the stable regime, each Gaussian 
peak has an amplitude which fits nicely to 



(T + To)" 



(46) 



where the fit parameters vary with the coupling constant 
R. In the i? — )■ limit, we know from Eq. (|40p that the 
first structure peak decays exponentially {Aq = 1, a = 0, 
E = TTr/6'r/S{Qi)) and that the prepeak is non-existant 
{Aq = 0). In the i? — ^ i?* limit, we find that both peaks 
tend toward pure power law decays (i? — >■ 0). The peak 
width can be nicely fit to 



B, = Bo{T + ToY 



(47) 



where /3 is approximately 1 below the critical coupling, 
but decreases very rapidly as R R* . The peak centers 
Qx drift slightly at early times, but quickly reach a fixed 
asymptotic value. Example fits for A and B can be seen 
in Fig. m 

In the unstable regime, peak amplitudes initially de- 
crease then increase without bound and the peak widths 
narrow. At early times and away from Tmin, the ampli- 
tudes may again be fit by Eq. (|46l) and the widths by 
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FIG. 7; Plots of the (a) amplitude and (b) width fit to Eqs. 
dm and (gZll for J7 = 0.55, R = 0.14. In (a) the fit is 
given by A{T) = Q.^Ze-^'^''^'^ /{T + Q.Q^mf -^^^ and in (b) 
by B{T) = 124.77T'' *°^ - 8.13. The data and fit are in- 
distinguishable except at the very earliest times. In (c), the 
structure peaks are shown (solid lines) with their Gaussian fits 
(dashed lines) at times T = 0.1, 0.25 and 0.5 (from largest to 
smallest amplitude). Note that the fit improves with time. 



One solution would be to extend the model by choos- 
ing an effective Hamiltonian known to support freezing. 
According to density functional theory, we can choose 
the effective Hamiltonian to be 



Xi 



(xi)ln(:^:i^)-<50(xi) 



-i J d'':Eid'^a;2<50(xi)Ci3(xi -X2)(50(x2) (48) 

which, after expanding the logarithm and rearranging, 
generates higher order terms compared to the Hamilto- 
nian used in this work: 

= ^y"d''xid'^X2(5(^(xi)x-l(xi -X2)<5</>(X2) 

+v J d^xi{6<j)ixi)f + u J d'^xi((5</)(xi))4 + ... (49) 
where 

(5(x) — n 



r'(x) = x-'(x) 



nl3 



(50) 



and the higher order coefficients are given by u = 
-1/6/3^2 and u = 1/I2(3n^. 

A full treatment of this model with the above static 
perturbations will be given elsewhere. 



D. Prepeak Behavior 



Eq. (j47p . However, the subsequent transition and growth 
pieces do not easily lend themselves to a universal fit, es- 
pecially at late times when the growth accelerates and 
the numerical solution breaks down. 

We find that these Gaussian peaks (in both the stable 
and unstable regimes) are tending toward delta functions 
as T — > oo; the peak widths are narrowing to zero and 
the amplitudes are time-dependent. This is the same 
behavior which was seen with the prepeak in RDMI. 



C. Unstable Regime 

The behavior of the model in the unstable regime is 
clearly unphysical. What is happening? At first, one 
might suspect that the system is undergoing a liquid to 
glass transition. However, as was shown in RDMI, an er- 
godic non-ergodic transition can be found for this model 
at lowest order, but is not sustained at higher order. It 
is therefore more likely that the system is instead simply 
freezing from a liquid to a solid. Since the model is for 
a system in equilibrium, as the density and coupling are 
increased, the system may wish to nucleate. Our model, 
however, is not equipped to handle such a transition and 
so the peaks grow and are not properly stabilized. 



We have seen that the prepeak is not always visible 
near and beyond the transition from stable to unstable. 
Let us explore the behavior and significance of the pre- 
peak in greater depth. 



1. Extraction of 



Information 



The prepeak is not always easily identifiable from the 
plots, but is hidden by other features in the data (e.g., 
the tails of the Q=0 peak and/or the first structure peak) 
and is of relatively small amplitude. To disentangle the 
prepeak from the other behavior, recall that we know the 
R = solution exactly (Eq. (|40l) ) and that the solution 
includes no prepeak and has a purely exponential decay 
for the structure factor peaks. 

Let us divide out the R = contribution to f{Q,T), 



HQ,T) 



f{Q,T) 
f{Q,T)\n^o' 



(51) 



Though f{Q,T) depends on R through an integro- 
differential equation and cannot be separated into purely 
R — and R ^ contributions in either an additive or 
multiplicative fashion, the method has the potential to 
yield some information since we have seen that the late 
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FIG. 8: (a) Plot of f(Q, T) for = 0.45, R = 0.42. (b) 
Logarithmic plot of h{Q,T) (given by Eq. ([ST}) in the pre- 
peak region at the same ry and R values. Times are T=0.02, 
0.05 and 0.1 (from smallest to largest peak amplitude). Note 
that while no prepeak is visible in f(Q,T), a rapidly growing 
Gaussian appears in h(Q,T) 



time behavior is a simple sum of Gaussian peaks, 
explicitly, we expect 



h{Q,T). 



foQ,T + fppiQ,T) + h{Q,T) 
/o(g,T) + /{(Q,r) 



More 



(52) 



where the prime simply designates that while the form 
for the peak is the same, the amplitudes and widths differ 

As a test of the above proposed methods, let us look 
at the plot of f{Q, T) for -q = 0.45, R = 0.42 (Fig. [8(a) I. 



We see that we are in the unstable region since the first 
structure peak is growing with time. However, we see 
no sign of any prepeak. Let us divide out the -q = 0.45, 
R = val ues from our original f{Q,T) to get h{Q,T) 
(Fig. |8(b)[ ). Now we find a very clear prepeak at Q ~ 3.40 
which grows very rapidly in amplitude and narrows in 
width with time along with a smaller peak in the location 
of the first structure peak which also grows and narrows. 

If we repeat this technique with other combinations of 
r] and R we find that regardless of the parameters, a pre- 



peak centered near Q « 3.40 can always be found. Also, 
whereas the prepeak and first structure peak centers ex- 
hibit some early time transient drift when the raw results 
are plotted, these new plots show fixed peak centers. 

Let us now try to fit these peaks to Gaussians. The 
peak amplitudes can be fit to 



A{T)^AoiT + Toy 



„ET 



(53) 



where a is positive for the first structure peak and neg- 
ative for the prepeak. The width is given by 

B{T) = B^{T + nf 



(54) 



where j3 is approximately 1 for both peaks. 



2. Connection to Experimental and Simulation Results 

While a prepeak developed in both the structureless 
approximation of RDMI and our own more realistic cal- 
culations here, there is little experimental evidence for 
the existence of such a peak in hard-sphere liquids. For 
molecular systems, a number of experiments reveal pre- 
peaks below the first structure peak [23l - |25| . but this 
work is likely unrelated to our model since the peaks 
can generally be explained through physical mechanisms 
unique to molecular systems (e.g. hydrogen bonding 
leading to clustering, molecular reconfigurations, etc.) 
and since these prepeaks generate static structure pre- 
peaks instead of the dynamic structure we see in our 
model. Work with molecular dynamics simulations on 
the other hand, has shown some dynamically generated 
prepeaks |26l - [29| . but again this work is on molecular 
systems and can be explained by uniquely molecular 
dyanamics, e.g. translation-rotation coupling. 

While no evidence for prepeaks in hard-sphere or 
monatomic systems has been reported, we show here that 
at higher densities the prepeak was not readily visible in 
the dynamic structure factor except in the (non-physical) 
long-time breakdown. It is therefore conceivable that a 
prepeak similar to that found here may in fact exist in 
experimental results, but has simply not been detected 
and extracted. 

How could a prepeak hidden in experimental data be 
extracted? In the previous section, we showed that hid- 
den prepeaks appear if one plots the raw data divided 
by the theoretically determined zero-coupling data as in 
Eq. (1511) . This form, however, is not useful to an experi- 
mentalist, so let us rewrite it in terms of density n, tem- 
perature j3 and the physical diffusion coefficient. Dp = D 
- parameters which would be known or could be deter- 
mined - and the static and dynamic structure factors 
{S{q) and C(g, t) respectively) - which would be mea- 
sured. We then have 



h{q,t) 



nS{q) 



exp 



q^Dpt \ 
nS{q)f3j- 



(55) 



For such a method to yield useful information, ex- 
tremely good data resolution for the dynamic structure 
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factor at small wavenumbers in the prepeak region would 
be needed, as would small uncertainties on the quantities 
going into Eq. (|55|) . Poor resolution and large uncertain- 
ties would easily wash out the small values that one is 
trying to extract. 

V. CONCLUSIONS 

We began with the random diffusion model introduced 
in RDMI which describes a system undergoing diffusive 
dynamics with a density-dependent bare diffusion coeffi- 
cient. While the model had previously been studied in 
the structureless approximation where the short distance 
degrees of freedom had been integrated out, we wished to 
study the model from a more realistic view that included 
such features. By relaxing the delta function constraint 
of the transport matrix and introducing the smoothing 
functions foil) and fi{q), we were able to introduce a re- 
alistic structure factor back into the problem while still 
keeping the wavenumber integrals finite. The resulting 
kinetic equation then depended on two variable param- 
eters, the dimensionless density rj and the coupling con- 
stant R. 

In our investigations, we came to the following conclu- 
sions: 

1. When the coupling constant is increased, there is a 
significant slowing down of the system; exponential 
decay gives way to algebraic decay. 

2. For large enough coupling there is a transition 
where the system goes from stable to unstable. All 
peaks grow without bound and the solution be- 
comes unphysical. The critical coupling required 
to reach this transition decreases as the density in- 
creases. 

3. It is possible this instability may be cut off by 
including terms in the model which allow for nu- 
cleation. We propose terms motivated by density 
functional theory and the results of this model ex- 
tension will be pursued in another work. 

4. Near but below the transition, a new peak, termed 
the prepeak, sometimes can be seen between the 
Q — peak and the first structure peak. When 
it can be seen, its amplitude is less than that of 
the first structure peak. Above the transition, the 
prepeak always appears, though it may only be seen 
at late times when the peak amplitudes are growing 
without bound. 

5. Near but below the transition, both the prepeak 
and first structure peak narrow to delta functions, 
but with algebraically decaying amplitude. 

6. It is possible to isolate the prepeak by dividing out 
the R = solution. In such a case we find the pre- 
peak and structure peaks to be well separated with 



growing amplitudes and narrowing widths. This 
technique can separate out the prepeak even when 
it cannot be seen in the raw data and we find that 
the prepeak is always centered near Q « 3.40 re- 
gardless of coupling or density. 

7. By the above technique, one should be able to 
take an experimentally determined dynamic struc- 
ture factor and isolate the prepeak using only ex- 
perimentally measured or controlled parameters. 
Such an isolation might reveal previously unde- 
tected prepeak information and give indirect mea- 
surements of parameters of the system not usually 
accessible such as the the coefficients of the terms 
in the bare diffusion coefficient. 

These results confirm that the features seen in RDMI 
were not due to the nature of the structureless approxi- 
mation, but were indeed inherent features of the model. 
In addition, the model is now sufficiently realistic to in- 
vite comparison with experimental results. While the 
prepeak seen in our model has not been reported exper- 
imentally, we find that at reasonably accessible param- 
eter values, it is often hidden; we provide a method by 
which a hidden prepeak might be extracted. The model, 
however, still lacks a sufficient mechanism to allow for 
freezing. Future work introducing static perturbations 
to the Hamiltonian may rectify this. 
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Appendix A: Method of Numerical Solution 

To solve the kinetic equation given by Eqs. ([50]) and 
(1311) . we use a simple Euler-type integration scheme to 
find f{Q,T) over a discrete set of wavenumbers Qi and 
times Tj . The Q values are equally spaced, but the T val- 
ues are chosen using an adaptive step routine and stored 
in an array. In various integration steps, values of f(Q,T) 
at points off this lattice are sometimes required. Since 
f(Q,T) varies smoothly, a simple interpolation scheme is 
used to estimate such values. Let us address the integra- 
tion method, the adaptive step routine and the interpo- 
lation scheme in turn. Methods were adapted from Ref. 

n. 

1. Integration Method 

The kinetic equation (Eq. (ISOjl Hs presented in the sim- 
ple form 

dfiQX) 

— — =g[Q,T) (Ai) 



12 



where 



giQ,T) 



'S{Q) Qrj 



f{Q,T) 



-R^ ( dSN{Q,T~ S)f{Q,S). (A2) 

"'0 



If we ignore the intricacies of ^(Q, T) for the moment, we 
can advance our solution to f{Q,T) one time step at a 
time (in the normal Euler fashion) as 

= /(g.,r,) 

+{T,+^-T,)g{Q^,T,) (A3) 

using the initial conditions that fiQi^T^) = 1. 

The function g{Q,T) includes a time integral over the 
memory kernel which can be decomposed into a similar 
differential equation, 



dhiQ,T,S) 
dS 



R^N{Q,T-S)f{Q,S), (A4) 



and can thus similarly be solved. Note now that since 
the time steps are not even, interpolation may be needed 
to find values of functions off the lattice. This will be 
addressed in detail below. 

At the next nested level, the memory function is itself 
an integral over wavenumber. In three dimensions, dK = 
K^dKdsmOdcj). Integrating over cj) and expressing dot 
products as |Q • K| = y/Q'^ + K'^ — 2QKcos9 reduces 
the integral to two dimensions, which are again solved 
by the Euler method using interpolations. 

Thus, to summarize, as we advance f{Q,T) from one 
time step to the next, we must first advance the memory 
function array by one time step, then perform the time 
integration over it to get the function g{Q, T) and finally 
use g{Q,T) to fully advance f{Q,T). 



1. Adaptive Step Routine 

We desire an adaptive time step size since we expect 
our solutions to sometimes offer smooth decays to zero 
and other times offer rapid growth and eventual instabil- 
ity. Let us describe the routine. 

When choosing step sizes, there is always a trade- 
off between maximizing computation speed (larger step 
sizes) and minimizing error (smaller step sizes). To quan- 
tify this trade off, let us imagine advancing a function 
y{x) from initial value y{xo) up to y{xf) in two differ- 
ent ways. In the first method, let us take two steps of 
size hi = {x f — xo)/2 and call the result Yi. In the sec- 
ond method, let us instead take one large step such that 
/i2 — 2hi and call that result ¥2- The difference between 
these two results is then AY = \Yi — Y2\. If one now 
specifies a desired difference between these values AFq 



- which in some sense quantifies the trade off between 
accuracy and speed described above - then we can infer 
that the best choice for the next step is given by 



^new — ^1 



AYo 
AY' 



(A5) 



We see that if the actual difference is smaller than our 
desired difference, we can afford to increase the step size, 
whereas if the actual difference is larger, we decrease our 
step size. 

In our numerical integration, we specify an initial 
time step and tolerance AYq. After each time step, 
the program determines where in Q-space the fastest 
growth/decay is occurring (where one assumes AY will 
be the largest) and then computes the two test values 
described above at that point; the program integrates up 
from f{Q,Tj-i) to f{Q,Tj+i) first with two steps of the 
current size and then again with one step at twice that. 
From these points, the next step is determined. 

Two limiting actions are taken to prevent run-away 
changes from occurring. First, a safety factor of 0.9 is 
added to the right hand side of Eq. (IA5|) . This serves to 
underestimate the new step size so that the adjustments 
always err on the side of smaller error rather than faster 
speed. Second, the new step size is never allowed to 
be more than four times the previous step size. This 
prevents dramatic changes in step size. 

The choices of these limiting values (0.9 and four times) 
as well as our choice for the tolerance AYq = 10~^ were 
chosen through rough trial and error. 



3. Interpolation Algorithm 

Integrations over the memory kernel require sampling 
the function at times not on the previously computed 
Qi, Tj lattice. In such cases, a two-dimensional bilinear 
interpolation was used. If one wants to know the value 
of f{Q,T) and knows the value at surrounding points 
(Qi,Ti), (Qi,T2), iQ2,Ti) and (Q2, ^2) where Qi <Q< 
Q2 and Ti < T < T2, then we can estimate the value of 
f{Q,T) as 



/(Q,T) = (i-i)(i-«)/(Qi,ri) + t(i-u)/(Q2,ri) 

+tuf{q2,T2) + {l-t)uf{Qi,T2) (A6) 



where we define 



and 



i= (Q-Qi)/(Q2-Qi) 



u={T-Ti)/iT2~Ti). 



(A7) 



(A8) 
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